package com.chaintor.demo.math;

import org.apache.commons.math3.analysis.differentiation.DerivativeStructure;
import org.apache.commons.math3.analysis.interpolation.SplineInterpolator;
import org.apache.commons.math3.analysis.polynomials.PolynomialFunction;
import org.apache.commons.math3.analysis.polynomials.PolynomialSplineFunction;
import org.apache.commons.math3.distribution.NormalDistribution;
import org.apache.commons.math3.util.FastMath;

/**
 * deltaSimple的三次样条曲线函数
 */
public class DeltaSimpleFunction {
    private NormalDistribution normal = new NormalDistribution();
    private double atm,bf10,bf25,rr10,rr25;
    private double d10,d25,d50,d75,d90;
    private double ds10,ds25,ds50,ds75,ds90;
    private double[] knots;
    private PolynomialSplineFunction dspsfunc;       //三次样条函数
    private double[][] pc;
    private PolynomialSplineFunction derivativeFunc; //三次样条函数的导函数

    /**
     * 使用ATM报价构造
     * @param datm
     * @param drr25
     * @param dbf25
     * @param drr10
     * @param dbf10
     * @param te_to_t
     */
    public DeltaSimpleFunction(double datm,double drr25,double dbf25,double drr10,double dbf10,double te_to_t){
        //接收ATM报价
        this.atm = datm/100.0;
        this.bf10 = dbf10/100.0;
        this.bf25 = dbf25/100.0;
        this.rr10 = drr10/100.0;
        this.rr25 = drr25/100.0;
        //ATM报价到delta关键点
        this.d50 = this.atm;
        this.d10 = this.atm + this.bf10+0.5*this.rr10;
        this.d25 = this.atm + this.bf25+0.5*this.rr25;
        this.d75 = this.atm + this.bf25-0.5*this.rr25;
        this.d90 = this.atm + this.bf10-0.5*this.rr10;
        //deltaSimple关键点
        this.ds50 = normal.cumulativeProbability(-(d50* FastMath.sqrt(te_to_t/365.0)/2.0));
        this.ds10 = normal.cumulativeProbability(normal.inverseCumulativeProbability(0.1)-(d10*FastMath.sqrt(te_to_t/365.0)/2.0));
        this.ds25 = normal.cumulativeProbability(normal.inverseCumulativeProbability(0.25)-(d25*FastMath.sqrt(te_to_t/365.0)/2.0));
        this.ds75 = normal.cumulativeProbability(normal.inverseCumulativeProbability(0.75)-(d75*FastMath.sqrt(te_to_t/365.0)/2.0));
        this.ds90 = normal.cumulativeProbability(normal.inverseCumulativeProbability(0.90)-(d90*FastMath.sqrt(te_to_t/365.0)/2.0));

        SplineInterpolator spline = new SplineInterpolator();
        double x[] = {ds10,ds25,ds50,ds75,ds90};
        double y[] = {d10,d25,d50,d75,d90};
        this.knots = x;
        dspsfunc = spline.interpolate(x,y);
        derivativeFunc = dspsfunc.polynomialSplineDerivative();

        pc = getC();

    }

    /**
     * 使用参数构造
     * @param knots
     * @param c
     */
    public DeltaSimpleFunction(double[] knots,double[][] c){
        PolynomialFunction ps[] = new PolynomialFunction[knots.length-1];
        for( int i=0;i<knots.length-1;i++){
            ps[i] = new PolynomialFunction(c[i]);
        }

        this.knots = knots;
        this.pc = c;
        dspsfunc = new PolynomialSplineFunction(knots,ps);
        derivativeFunc = dspsfunc.polynomialSplineDerivative();
    }

    /**
     * 按delte调整好的报价构造
     * @param type 为区别ATM报价添加的参数，随便填
     * @param datm
     * @param dd10
     * @param dd25
     * @param dd75
     * @param dd90
     * @param te_to_t
     */
    public DeltaSimpleFunction(String type,double datm,double dd10,double dd25,double dd75,double dd90,double te_to_t){
        this.atm = datm/100.0;

        this.d50 = datm/100.0;
        this.d10 = dd10/100.0;
        this.d25 = dd25/100.0;
        this.d75 = dd75/100.0;
        this.d90 = dd90/100.0;

        //deltaSimple关键点
        this.ds50 = normal.cumulativeProbability(-(d50* FastMath.sqrt(te_to_t/365.0)/2.0));
        this.ds10 = normal.cumulativeProbability(normal.inverseCumulativeProbability(0.1)-(d10*FastMath.sqrt(te_to_t/365.0)/2.0));
        this.ds25 = normal.cumulativeProbability(normal.inverseCumulativeProbability(0.25)-(d25*FastMath.sqrt(te_to_t/365.0)/2.0));
        this.ds75 = normal.cumulativeProbability(normal.inverseCumulativeProbability(0.75)-(d75*FastMath.sqrt(te_to_t/365.0)/2.0));
        this.ds90 = normal.cumulativeProbability(normal.inverseCumulativeProbability(0.90)-(d90*FastMath.sqrt(te_to_t/365.0)/2.0));

        SplineInterpolator spline = new SplineInterpolator();
        double x[] = {ds10,ds25,ds50,ds75,ds90};
        double y[] = {d10,d25,d50,d75,d90};
        this.knots = x;
        dspsfunc = spline.interpolate(x,y);
        derivativeFunc = dspsfunc.polynomialSplineDerivative();
        pc = getC();

    }

    /**
     * 计算
     * @param x
     * @return
     */
    public double value(double x){
        if (x<knots[0]){
            return pc[0][0]+pc[0][1]*(x-knots[0])+pc[0][2]*(x-knots[0])*(x-knots[0])+pc[0][3]*(x-knots[0])*(x-knots[0])*(x-knots[0]);
        }
        if (x>knots[4]){
            return pc[3][0]+pc[3][1]*(x-knots[3])+pc[3][2]*(x-knots[3])*(x-knots[3])+pc[3][3]*(x-knots[3])*(x-knots[3])*(x-knots[3]);
        }
        return dspsfunc.value(x);
    }

    /**
     * 导数计算
     * @param x
     * @return
     */
    public double dvalue(double x){
        if (x<knots[0]){
            return pc[0][1]+2.0*pc[0][2]*(x-knots[0])+3.0*pc[0][3]*(x-knots[0])*(x-knots[0]);
        }
        if (x>knots[4]){
            return pc[3][1]+2.0*pc[3][2]*(x-knots[3])+3.0*pc[3][3]*(x-knots[3])*(x-knots[3]);
        }
        return derivativeFunc.value(x);
    }


    /**
     * 返回各段的方程系数
     * @return
     */
    public double[][] getC(){
        PolynomialFunction[] ps = dspsfunc.getPolynomials();
        double[][] c = new double[ps.length][];
        for(int i=0;i<ps.length;i++){
            c[i] = ps[i].getCoefficients();
        }

        return c;
    }

    /**
     * 返回各段的端点值
     * @return
     */
    public double[] getKnots(){
        return dspsfunc.getKnots();
    }
}
